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Abstract. In this paper we give an introduction to the numerical density matrix 
renormalization group (DMRG) algorithm, from the perspective of the more general 
matrix product state (MPS) formulation. We cover in detail the differences between 
the original DMRG formulation and the MPS approach, demonstrating the additional 
flexibility that arises from constructing both the wavefunction and the Hamiltonian in 
MPS form. We also show how to make use of global symmetries, for both the Abclian 
and non- Abclian cases. 
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1. Introduction 

The DMRG algorithm was introduced by Steven White [TJ, as an algorithm for 
calculating ground state properties of principally one-dimensional strongly correlated 
systems in condensed matter physics. The connection between DMRG and matrix 
product states [21 El HI [5] (also known as finitely correlated states) was first made 
by Rommer and Ostlund [6J, who identified the thermodynamic limit of DMRG 
with a position-independent matrix product wavefunction. Although DMRG had 
already proven itself to be useful empirically, this was an important step in rigorously 
establishing the physical basis of the algorithm due to the concrete and easy to 
manipulate form of matrix product states. Work on the spectra of density matrices 
[7] , later formulated as scaling of the von Neumann entropy [8] has placed the algorithm 
on a firm footing, showing that the required computational effort (realized via the basis 
dimension m) is essentially a function of the entanglement of the wavefunction [9] , which 
for one- dimensional ground-states scales at worst logarithmically with the system size 

ma. 

Computationally, MPS algorithms came to the fore with the assistance of a quantum 
information perspective, leading to algorithms for periodic boundary conditions [TT], and 
finite temperature algorithms based on density operators [131 EEl EE]- At around the 
same time, methods for simulation of real time evolution were developed in DMRG 
[T5l [TB]. which can also benefit from MPS formulations [TT] . 

The common theme of MPS approaches is to allow algorithms that operate on 
multiple, distinct wavefunctions at the same time. This is possible in the original 
formulation of DMRG only by constructing a mixed effective Hilbert space that is 
weighted appropriately to represent all of the relevant states simultaneously. This is 
inefficient, as the algorithms typically scale as 0(m 3 ) (or up to 0(m 5 ) for periodic 
boundary conditions [TT]) in the number of basis states m, so increasing m so as to 
represent multiple states in the same basis is typically much slower than performing 
separate operations on each basis. In addition, the mixed basis approach lacks flexibility. 
While traditional DMRG programs calculate the wavefunction and a few (often 
predetermined) expectation values or correlation functions, if instead the wavefunction 
is calculated in the MPS representation of Eq. (0Q) it can be saved for later use as an input 
for many purposes. Perhaps the simplest such operation beyond the scope of traditional 
DMRG is to calculate the fidelity, or overlap between the wavefunctions obtained from 
separate calculations. In the MPS formulation, this calculation is straightforward. 
Nevertheless the determination of the scaling function for the fidelity of finite-size 
wavefunctions for different interaction strengths, provides a new tool for investigating 
phase transitions and crossover phenomena [1_H HH [20] . Indeed, due to the simplicity 
of the calculation the fidelity is likely in the coming years to be the first choice for 
quantitatively determining critical points. Similar measures of entanglement, such as 
the concurrence and single- and two-site entropy [2TJ, [22], are also straightforward to 
calculate, hence the MPS formalism allows us to apply directly the emerging tools of 
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quantum information to the study of realistic systems in condensed matter physics. An 
alternative measure, the Loschmidt Echo [23J is important because, unlike many of the 
quantum information theoretic measures, this is directly accessible in experiments while 
showing the rich behavior of the simpler measures. The Loschmidt Echo is more time- 
consuming to measure numerically as it requires a full time evolution simulation rather 
than a direct measurement, nevertheless it is well within the current state of the art 

In this paper, we focus on the case of open boundary condition matrix product 
states. This does not preclude calculation of periodic systems, however the entanglement 
of such periodic states is increased such that in the large L limit (where L is the lattice 
size), the number of states kept tends to the square of that required for open boundary 
conditions [25J. Algorithms exist for periodic boundary conditions [TT] and infinite 
systems [26] (not to be confused with the 'infinite-size' DMRG algorithm), and the 
basic formulas introduced here carry over to these cases, but we do not describe the 
specific algorithms here. In Sec. [2} we introduce the basic formulation of matrix product 
states, and formulas for the fidelity. Sec. [3] is devoted to a new approach, whereby we 
construct the Hamiltonian operator itself as an MPS, with many advantages. We cover 
some remaining details of the DMRG algorithm in Sec. HI before discussing in detail 
the use of Abelian and non-Abelian quantum numbers in Sec. El We finish with a few 
concluding remarks in Sec. El including some observations on finite temperature states. 

2. Matrix Product States 

We denote an MPS on an L-site lattice by the form 

Tr^A Sl A S2 ---A^ | |a 2 )®...® \s L ) , (1) 
{*} 

The local index Sj represents an element of a local Hilbert space at site i. The two 
important cases we cover here are when Sj runs over a <i-dimensional local basis for 
a wavefunction |sj), in which case we refer to this as a matrix product wavefunction 
(MPW), or Si is a (i 2 -dimensional local basis for all operators acting on a local site, 
which we refer to as a matrix product operator (MPO). In this paper, we use MPS 
for a generic state irrespective of the form of the local space, and use MPW or MPO 
as necessary when the distinction between wavefunctions and operators is important. 
In general, the matrix product form can also represent periodic [11] and infinite (non- 
periodic) states [26], but here we use only the open-boundary form equivalent to the 
wavefunction obtained by DMRG [6J. To enforce this boundary condition, we require 
the left-most matrix A Sl to be 1 x m dimensional, and the right-most matrix A SL to be 
m x 1. Here we have introduced m as the basis size, or dimension of the matrix basis 
of the v4-matrices. This quantity is often denoted D, or sometimes Xi 111 the quantum 
information literature, but we emphasize it is exactly the same quantity in all cases. In 
general m is position dependent, as we do not require the ^-matrices to be square even 
away from the boundary. Because of the 1-dimensional basis at the boundaries we can 
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regard the MPS wavefunction to be a sequence of operators attached to left and right 
(or outgoing and incoming) vacuum states. This makes the operator product in Eq. ([T|) 
an ordinary number, so the trace operation can be dropped. 

2.1. Orthogonality constraints 

In practice, a MPS state with no particular constraints on the form of the A-matrices 
is numerically difficult to handle. We are always free to insert some product of a non- 
singular m x m operator X and its inverse X~ l in the middle of our MPS, thus we 
can apply an arbitrary transformation to the matrix basis of an A-matrix, as long as 
we make the corresponding transformation to its neighbor. Using this freedom, we can 
transform the A-matrices into a form where they are orthonormalized, that is, we prefer 
that they satisfy one of two possible constraints, 

E s A s A st = 1 (right-handed) 
J2 S A S ^A S = 1 (left-handed) 

States satisfying these conditions are orthonormalized in the sense that if all A-matrices 
to the left of some matrix A Sn are orthogonalized in the left-handed sense, then the 
basis on the left-hand side of A Sn is ortho normal [ie the identity operator in the 
effective Hilbert space is trivial). Conversely, if all A-matrices to the right of A Sn are 
orthogonalized in the right-hand sense, then the basis on the right-hand side of A Sn is 
orthogonal. Usually, we want both these conditions to be true simultaneously. Note 
that it is not, in general, possible for all of the A-matrices (including A Sn itself) to be in 
orthonormal form at the same time. There are several ways of transforming an arbitrary 
MPS into this normalized form. Two ways that we consider here are the singular value 
decomposition (SVD), and the related reduced density matrix, as used in DMRG pQ. 
The simplest, and in principle the fastest, is the SVD, well-known from linear algebra 
[27]. For example, for the left-handed orthogonality constraint on Afj, where we have 
re-inserted the matrix indices i,j, we consider s,i to be a single index of dimension 
dm, giving an ordinary dm x m dimensional matrix, and carry out the singular value 
decomposition, 

4 = EV«(^)«. ( 3 ) 

kl 

where U is column-orthogonal, U'U = 1, and is row-orthogonal, V*V — 1. D 
is a non-negative diagonal matrix containing the singular values. This form coincides 
with the Schmidt decomposition, where D gives the coefficients of the wavefunction in 
the Schmidt basis [21 J. The matrix U therefore satisfies the left-handed orthogonality 
constraint, so we use this as the updated A-matrix, and multiply the A-matrix on the 
right by DV\ This implies that the ^4-matrix on the right is no longer orthonormalized 
(even if it was originally), but we can apply this procedure iteratively, to shift the non- 
orthogonal A-matrix to the boundary - or even beyond it - at which point the lxl 
A-matrix coincides with the norm of the wavefunction. An important point here is 
that we can choose to discard some of the states, typically those that have the smallest 
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singular value. This reduces the matrix dimension m, at the expense of introducing 
an approximation to our wavefunction, such that the squared norm of the difference 
of our approximate and exact wavefunctions is equal to the sum of the squares of the 
discarded singular values. Note however that the singular values only correspond to 
the coefficients of the Schmidt decomposition if all of the remaining A-matrices are 
orthogonalized according to Eq. (j2J). If this is not the case, the singular values are not 
useful for determining which states can be safely discarded. 

Alternatively, we can construct the reduced density matrix, obtained by tracing 
over half of the system. This is achieved by 

Pif = X! A S ik*A S j k , (4) 

k 

which is a dm x dm matrix, with m eigenvalues coinciding with the values on the 
diagonal of D 2 , and the remaining eigenvalues are zero. Again, the eigenvalues are only 
meaningful if the remaining A-matrices are appropriately orthogonalized. The utility of 
the density matrix approach over the SVD, is that we can introduce mixing terms into 
the density matrix which can have the effect of stabilizing the algorithm and accelerating 
the convergence, which is further discussed in Sec. HJ 

The overlap of two MPS is an operation that appears in many contexts. For 
wavefunctions this gives the fidelity of the two states, and for operators this is equivalent 
to the operator inner product (A\B) = TrA^B which induces the Frobenius norm. 
Direct expansion of the MPS form yields, 

(A\B) = ]T(Tr A S1 *A S2 * ...){Ty B Sl B S2 ...) 
= Tr J2( AS1 * ® B ' S1 )( AS2 * ® B ' S2 ) ■ ■ ■ ■ 

{Si} 

Due to the open boundary conditions, the direct product E\ = A Sl * ® B Sl reduces to 
an ordinary m x m matrix, after which can construct successive terms recursively, via 

E n = J2 Atf K-i B£ , (6) 

8 n 

with an analogous formula if we wish to start at the right hand side of the system and 
iterate towards the left boundary, 

F n = £ A£ F n+1 B^ . (7) 

For the purposes of numerical stability, it is advisable for the A- and S-matrices to be 
orthogonalized in the same pattern, that is, E'-matrices are associated exclusively with 
the left-hand orthogonality constraint and F-matrices are associated with the right-hand 
orthogonality constraint. If we iterate all the way to the boundary, the E- (or F-) matrix 
ends up as a 1 x 1 matrix that contains the final value of the fidelity. Alternatively, we 
can iterate from both ends towards the middle and calculate the fidelity as TrEFK 
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2.2. Local updates 

The key notion of the matrix product scheme is that of local updates] that is, we 
modify, typically though some iterative optimization scheme, one (or perhaps a few) A- 
matrix while keeping the remainder fixed. A useful and flexible alternative is the center 
matrix formulation, where, instead of modifying an A-matrix directly, we introduce an 
additional matrix C into the MPS, 

A* A" ■ ■ ■ A Sn C A Sn+1 ■ ■ ■ A SL . (8) 

This allows us to preserve orthogonality of the matrices at all times; matrices A Si for 
i < n are normalized always according to the left-handed constraint, and matrices for 
i > n are normalized according to the right-handed constraint. We directly modify only 
the matrix C which simplifies the local optimization problem as C is just an ordinary 
matrix. To introduce the local degrees of freedom, say for the | s n ) states, we expand 
the basis for A Sn . That is, we replace the mxm dimensional matrices A Sn with m x dm 
matrices A' Sn , given by 

A!* = 5 id+Snd , (9) 

and introduce the dm x m dimensional center matrix 

= AS , (io) 

with j = id + s n running over dm states. This doesn't change the physical wavefunction, 
as A' Sn C = A Sn . Similarly, we can expand the basis for the A-matrix on the right side of 
C, to give the effect of modifying either a single A-matrix, or two (or more) at once. In 
the center matrix formulation, the singular value decomposition required for truncating 
the basis is simply the ordinary SVD on the matrix C = UDV\ and we multiply (for a 
right-moving iteration) A Sn U, which preserves the left-handed orthogonality constraint, 
and DV^A Sn+1 which is not orthogonal, but becomes so when we again expand the basis 
to construct the new C matrix. The density matrix in the center matrix formulation 
is simply p = CC^ or p = C^C for left and right moving iterations respectively. For 
readers already familiar with DMRG, the center matrix corresponds exactly with the 
matrix form of the superblock wavefunction |29j . 



2. 3. Matrix product arithmetic 

The utility of the MPS approach is realized immediately upon attempting manipulations 
on the wavefunction Eq. (Q]). Suppose we have two distinct MPS, defined over the same 
local Hilbert spaces, 

\V A ) = A^A S ^...A S ^ |si)|s 2 >... 

|* fl ) = B Sl B s ' 2 . . . B SL |si)|s 2 )... 1 ' 

The superposition \ipc) — \^a) + \ ^b) is formed by taking the sum of the matrix 
products, A Sl A S2 . . . + B 81 B S2 . . ., which can be factorized into a new MPS, \^c) — 
C Sl C S2 with 

C si = a 3 * © B Si . (12) 
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To preserve the one-dimensional basis at the boundaries, the direct sum is replaced by 
a concatenation of columns or rows, for the left and right boundary respectively. This 
procedure results in an MPS of increased dimension, mc = rriA + tub- Thus, after 
constructing these matrices we need to re-orthogonalize the state, and then we can, if 
necessary, truncate the basis size to a given truncation error, which is well defined here 
and measures the exact difference between the original and truncated wavef unctions. 
Alternatively, the normalized and truncated MPS | ^c) can be constructed directly, 
by calculating the overlap matrices E between \ ^a) and \ ^b)- From the .E-matrices 
introduced in Eq. (Q, we can construct directly the orthogonalized reduced density 
matrices of | \I/c) an d truncate the basis as required, in a single step. This approach has 
better computational scaling than the two-step procedure of first orthogonalizing and 
then truncating, especially when the number of MPS in the superposition is large. But 
in general, iterative optimization approaches, where we use a DMRG-like algorithm to 
optimize the overlap (9c\ (I^a) + \ ^b)), have even better performance scaling with 
large m or more states in the superposition. 

3. Operators 

A useful generalization of the MPS structure Eq. ([1]) is to use it to represent an operator 
(an MPO) instead of a wavefunction. This has been used before for calculating finite- 
temperature density matrices [12], but here we instead want to use this structure to 
represent the Hamiltonian operator itself. All Hamiltonian operators with finite-range 
interactions have an exact MPS representation with a relatively small matrix dimension 
M. For example, the Ising model in a transverse field has a dimension M = 3, and the 
fermionic Hubbard model has dimension M = 6. We use the capital letter to distinguish 
from the dimension of the wavefunction, m. Similarly, the local dimension of the upper 
index is denoted here by D, which is usually just equal to d 2 , but slightly complicated 
in the case of non-Abelian quantum numbers (see Sec. 15.21) . 
We denote an MPO by the form 

^ M S 'AM^--M S ^ |si)<si|® |4)(s 2 |--- , (13) 

where again we require that the first and last dimensions are M — 1, for open boundary 
conditions. 

The orthogonality constraint used previously for the MPS, Eq. ([2]), is not 
appropriate for Hamiltonian operators. When applied to an operator, the 
usual orthogonality constraints utilize the (Frobenius) operator norm, which scales 
exponentially with the dimension of the Hilbert space. With this normalization, 
components of an MPO Hamiltonian, such as the identity operator or some interaction 
term, tend to differ in magnitude by some factor that increases exponentially with the 
lattice size. Arithmetic manipulations on such quantities is a recipe for catastrophic 
cancellation [27] resulting in loss of precision. Mixing operators with a unitary 
transformation (for example 0[ = OiCos9 + C^sin^, Of, = —Oism9 + C^cos^), 
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will lead to a disaster if 0\ and O2 differ by a sufficiently large order of magnitude, 
~ fO 16 for typical double-precision floating point arithmetic. But such rotations are 
inevitable in the orthogonalization procedure because in general the operator inner 
product 02) — TrO|0 2 will not be zero. Instead we completely avoid mixing 
different rows/columns of the operator M-matrices, only collapsing a row or column if 
it is exactly parallel with another row or column. In this case, the actual norm of each 
component of the A- matrices is irrelevant, as they are never mixed with each other (but 
see also the discussion of the single-site algorithm in Sec. H]). For physical Hamiltonian 
operators this remains essentially optimal, with the minimum possible matrix dimension 
M. The only operators for which this orthogonalization scheme does not produce an 
optimal representation are operators that have a form analogous to an AKLT [28] state 
where the local basis states of the S = 1 chain are replaced by local operators. The 
resulting operator contains an exponentially large number of iV-body interactions for 
all N — > 00. We know of no physical Hamiltonians of this form. 

Given a Hamiltonian as a sum of finite-range interactions, it is possible to construct 
the operator M-matrices such that they are entirely lower (or upper) triangular matrices, 
thus in principle we can 'normalize' the matrices via some kind of generalized LU or 
QR decomposition. In practice we don't need to do this, as the Hamiltonian can easily 
be constructed in lower-triangular form from the beginning. Imposing, again without 
loss of generality, that the top-left and bottom-right components of the operator M- 
matrices are equal to the identity operator /, we can construct the sum of 1-site local 
terms H = YU-^-i 

position-independent MPS, 



which we regard as a 2 x 2 matrix, the elements of which are d x d dimensional local 
operators. For nearest-neighbor terms, H = J2iXiY i+ i, we have 



with the obvious generalization to iV-body terms. The direct sum and direct product 
of lower triangular matrices is itself lower triangular, thus this form can be preserved 
throughout all operations. For open boundary conditions, the left (right) boundary 
1 x M (or Mxl) matrices are initialized to (0, ... , 0, I) and (1,0, ... , 0) T respectively. 

The principal advantage of formulating the Hamiltonian operator (and indeed, all 
operators needed in the calculation) in this way that it can be manipulated extremely 
easily, amounting to a symbolic computation on the operators. This is in contrast 
to the ad hoc approach used in past DMRG and MPS approaches where the block 
transformations required for each operator are encoded manually, with limited scope 
for arithmetic operations. In particular, the sum of operators is achieved via Eq. (TI2"]) . 
Products of operators are achieved by matrix direct product; given MPO's A and B, 




(14) 



M 



( I \ 
Y 



(15) 



\0 X I J 
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the product C = AB is given by the matrices 

C s ' s = Y J AS ' t ® B ts ■ (16) 
t 

An implication of this is that the square of an MPO has a matrix dimension of at most 
M 2 , which, since M is usually rather small, means that it is quite practical to calculate 
expectation values for higher-order moments, for example of the variance 

a 2 = (H 2 ) — (H) 2 = ((H-E) 2 ), (17) 

which has been mentioned previously [31] as it gives a rigorous lower bound on the 
energy (although with no guarantee that this corresponds to the ground-state). In 
practice this lower bound is too wide to be useful in all but the simplest cases, but of 
more interest is the property that the variance is, to first order, proportional to the 
squared norm of the difference between the exact and numerical wavef unctions, and 
therefore also proportional to the truncation error [33] (see Sec. HJ). Thus, this quantity 
gives a quantitative estimate of the goodness of the wavefunction even for situations 
where the truncation error is not available. For our numerical MPS algorithms the 
variance takes the role of the precision e in numerical analysis [27] via e ~ \fo 2 . 

Of a similar form to the product of two operators, the action of an operator M on 
a wavefunction | A) gives a wavefunction | B) with matrix elements, 

B s ' = M s ' s <g> A s . (18) 

s 

The MPO formulation also gives a natural form for the evaluation of expectation values, 
similarly to the fidelity of Eq. (jSJ) , 

(A\M\B) = J2Tr E%F^ , (19) 

a 

where the E- and F-matrices now have an additional index a that represents the terms 
in the MPO M, with a recursive definition 

where again we can either iterate all the way to a boundary, at which point the a index 
collapses down to one- dimensional and the E a or F a are lxl dimensional matrices 
containing the expectation value, or we can iterate from both boundaries and meet at 
the middle, where our expectation value is given by the scalar product Eq. ([191) . 

Incidentally, given that the identity operators occur in a fixed location in the 
operator M-matrix (ie. at the top-left and bottom-right of the M-matrix) this fixes 
the index a of the reduced Hamiltonian and identity operators for the left and right 
partitions of the system. That is, in the application of the Hamiltonian E a , F a matrices 
to the wavefunction we are guaranteed that the a = 1 component of E a <g> F a ^ 
corresponds precisely to H L eg) I R , and the a = M component corresponds to It ® H R . 
Thus, even after an arbitrary series of MPO computations we can still identify exactly 
which component of the E, F matrices corresponds to the block Hamiltonian. This is 




Figure 1. Polynomial expansion for the relaxation of the impurity magnetization 
in the SIAM, up to order t 12 . The x symbols denote the impurity magnetization 
calculated via adaptive time DMRG, the dashed line is a guide to the eye. Parameters 
of the calculation were (in units of bandwidth), T — 0.1, U — 0.2, h n — 0.1, edo = 0.05, 
on a log-discretized Wilson chain with A = 1.8. At time t = 0, the Hamiltonian was 
switched to hi = and e^i = —0.1. 



useful for eigensolver preconditioning schemes [30], for example to change to a basis 
where the block Hamiltonian is diagonal. 

As an example of the utility of this approach, Fig. [1] shows the time evolution of 
the magnetization of the impurity spin in the single impurity Anderson model (SIAM), 
where the ground-state is obtained with a small magnetic field which is then turned 
off at time t = 0. The MPS operator approach readily allows the evaluation of the 
commutators required for a small t expansion of the expectation value of an observable 
in the Heisenberg picture, 

Mt) =A+^[H,A]- [If, [ff, A]] - ^[H, [H, [H, A]]] + ■■■ . (21) 

Since the number of terms in the repeated commutator will, in general, increase 
exponentially the accessible time-scales from such an expansion are clearly limited. 
Nevertheless this is a very fast and accurate way to obtain short-time-scale dynamics, 
and in this example 12 th order is easily enough to determine the T\ relaxation rate. 
For this calculation, the terms up to t 8 took a few seconds to calculate, while the t 10 
term took 6 minutes and the t 12 term took just over an hour, on a single processor 
2GHz Athlon64. This time was divided between calculating the MPO matrix elements 
(the dimension of which was just over 2500 at the impurity site), and calculating the 
expectation value itself. 
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4. DMRG 

We now have all of the ingredients necessary to construct the DMRG algorithm for 
determining the ground-state. Indeed, given the previous formulations, the DMRG 
itself is rather simple; using the center matrix formulation, we iteratively improve the 
wavefunction locally by using the center matrix C as input to an eigensolver for the 
ground-state of the Hamiltonian. The details of this calculation are precisely as for 
DMRG, already covered elsewhere [29 J. 

An important component of DMRG, which has been neglected in some matrix 
product approaches, is the truncation error. If only a single site is modified at a time, 
the maximum number of non-zero singular values is bounded by the matrix dimension 
m, thus the matrix dimension cannot be incrementally increased as the calculation 
progresses. Some way of avoiding this limitation is practically essential for a robust 
algorithm. The original DMRG formulation [1] solved this problem by modifying two 
A-matrices simultaneously, equivalent to expanding the matrix dimension for both the 
left and right A-matrices in Eq. (|Sj) so that the center matrix has dimension md x md. 
A scheme for single-site algorithms that avoids the limit on the singular values was 
introduced by Steven White [M], which uses a mixed density matrix constructed by 
a weighted sum of the usual reduced density matrix and a perturbed density matrix 
formed by applying the ^-matrices (on a right-moving sweep) or F a -matrices (on a 
left-moving sweep) of the Hamiltonian, 

P ' = p + cY,E a pE a ^ (22) 

a 

where c is some small factor that fixes the magnitude of the fluctuations. This solves 
nicely the problem of the bound on the number of singular values and introduces 
fluctuations into the density matrix that give the algorithm good convergence properties, 
often better than the two-site algorithm. A minor problem is that the scaling of the E a 
matrices is not well defined, in that we can scale the E^-matrices by an arbitrary M x M 
non-singular matrix X a i a , while at the same time scaling the F-matrices by X^ 1 . For 
one- and two-site terms, there is an 'obvious' scaling factor to use, whereby the scaling 
factors are chosen such that the (Frobenius) operator norm of the E and F-matrices 
are identical, but I don't know how this would apply more generally. An alternative 
that appears interesting is to apply the full Hamiltonian to a density operator for the 
full system, p' = p + cTr R H(p® I)H, constructed from the left (right) reduced density 
matrix and the right (left) identity operator, followed by a trace over the right (left) 
partition. In MPS form, this operation is 

f/ = p + c y £E a pE^G aP , (23) 

a,/3 

where G a p = TtF^F 13 is an M x M coefficient matrix. However, this scheme often 
fails; incorporating the G a p matrix reduces the fluctuations such that E a pE^ G a p 
differs little from p itself and the algorithm typically fails to reach the ground-state. 
The single-site algorithm [31] corresponds to choosing G a p = b~ a p. A useful compromise 
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appears to be using only the diagonal elements, such that G a p = 5 a/ 3 Tr F^F 13 , but this 
is surely not the last word on this approach. Both the two-site and mixed single-site 
algorithm inevitably result in a reduction in the norm of the wavefunction by truncating 
the smallest non-zero eigenvalues of the density matrix. The sum of the discarded 
eigenvalues, summed over all iterations in one sweep, is equal to the truncation error r, 
familiar from DMRG [T] (but note that it is common in the literature to quote an average 
or maximum truncation error per site). This quantity is useful in giving an estimate 
of the error in the wavefunction in this is, for a properly converged wavefunction, 
proportional to the norm of the difference between the exact ground-state and the 
numerical approximation. The presence of the truncation error explains why the bare 
single-site algorithm, despite having slightly better variational wavefunction than the 
two-site or mixed single-site algorithms [35J, converges much slower; the single site 
algorithm is a highly constrained optimization within an m-dimensional basis, whereas 
the two-site and mixed single-site algorithms are selecting the optimal m basis states out 
of a pool of a much larger set of states, namely the discarded states at each iteration 
(total ~ Lm states). While the notion of truncation error remains useful in MPS 
algorithms, for the purposes of error analysis we much prefer the variance Eq. ( TTTj) as 
being a direct measure of the accuracy of the wavefunction, independent of the details 
of a particular algorithm [33] . 

Low-lying excited states can be constructed using this algorithm. This has been 
done in the past in DMRG by targeting multiple eigenstates in the density matrix, 
but the MPS formulation allows a substantial improvement. Namely, it is easy 
to incorporate into the eigensolver a constraint that the obtained wavefunction is 
orthogonal to an arbitrary set of predetermined MPS's. That is, after constructing the 
MPS approximation to the ground-state, we can, as a separate calculation, determine the 
first excited state by running the algorithm again with the constraint that our obtained 
wavefunction is orthogonal to the ground-state. This is achieved by constructing 
the .E-matrices that project the set of states to orthogonalize against onto the local 
Hilbert space. These matrices are precisely those used in constructing the fidelity, Eq. 
([5]), thus given the center matrix of some state Cx, we project this onto the current 
Hilbert space C' x = ECxF^, and as part of the eigensolver, orthogonalize our center 
matrix against this state. This is a very fast operation, much faster than even a 
single Hamiltonian-wavefunction multiplication. So it is quite practical to orthogonalize 
against a rather large number of states, the practical limit is rather on numerical 
limitations in orthogonalizing the Krylov subspace in the eigensolver. If this is combined 
with an eigensolver capable of converging to eigenvalues in the middle of the spectrum 
(say, the lowest eigenvalue larger than some bound E ), then we need only a small 
number of states to orthogonalize against, say half a dozen states immediately above E 
in energy. In our numerical tests it seems to be rather common to skip eigenvalues, which 
is why we cannot simply orthogonalize against a single state. With a larger number of 
states to orthogonalize against, skipping eigenvalues is less of a problem as we are likely 
to recover the missing eigenstate on a later calculation. Using this approach, quantities 
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such as the level spacing statistics can be determined for system sizes far beyond exact 
diagonalization [33J. 

5. Quantum Numbers 

An important feature of matrix product states is that they can easily be constrained 
by quantum numbers representing the global symmetries of the Hamiltonian, as long as 
the symmetry is not broken by the spatial geometry of the MPS. For example, internal 
rotational symmetries such as U(l) and 577(2) [36] can be maintained exactly, but 
for a real-space lattice we cannot utilize the momentum in the same way because the 
representation itself violates this symmetrj||]. To achieve this, we impose a symmetry 
constraint on the form of the ^-matrices, so that they are irreducible tensor operators. 
That is, under a symmetry rotation the matrix A s for each local degree of freedom s 
transforms according to an irreducible representation D(j s ) of the global symmetry 
group. This is a very general procedure, that is applicable to essentially all MPS 
approaches and generalizations thereof. 

5.1. Abelian symmetries 

For Abelian symmetries, the representations are one-dimensional therefore the set of 
quantum numbers labeling the irreducible representations also forms a group, which we 
can write as 

D{j) ® D(k) = D(j + k) , (24) 

for two representations D(j) and D(k), where j + k denotes the group operation. Thus 
to incorporate Abelian symmetries into the algorithm we simply attach a quantum 
number to all of the labels appearing in the MPS, with the constraint that each A- 
matrix transforms irreducibly, so that the only non-zero matrix elements are 

A\, q q' = q + k, (25) 

where k, q', q are the quantum numbers attached to the local basis state and left and 
right matrix basis states respectively. We have suppressed here indices not associated 
with a quantum number, a convention which will be followed for the remainder of the 
paper. 

By convention, for our open boundary condition MPS we choose the right hand 
vacuum state to have quantum number zero. The symmetry constraint Eq. ( 1251) then 
implies that the quantum number at the left hand vacuum will denote how the state 
as a whole transforms (the target state, in DMRG terminology). This is the only real 
difference between DMRG and MPS algorithms, in that the DMRG convention is to 
construct both blocks starting from a scalar (quantum number zero) vacuum state, so 
that the superblock wavefunction is a tensor product of two ket states, 

\^f) = ip uv \ u) ® | v) , u + v = target , (26) 

uv 

I See also Ref. [37] for a real-space approach to constructing momentum eigenstates. 
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whereas for the MPS formulation the superblock wavefunction is represented by a scalar 
operator with a tensor product basis | u) ® (v\ with quantum numbers u = v. This means 
that, in contrast to the usual formulation of DMRG, the target quantum number is not 
encoded in the superblock but rather in the left vacuum state. A consequence of this 
is that DMRG is capable of representing simultaneously states with different quantum 
numbers, but an MPS is not. This is an important detail, for example in the calculation 
of dynamical correlations, as both the correction vector [38] and the similar DDMRG 
[39J algorithm require a basis optimized for both the ground-state | and the so-called 
Lanczos- vector where A is some operator that may not be scalar. However, the 

MPS formulation allows significant optimizations to these algorithms whereby the the 
calculation of the ground-state is decoupled from that of the Lanczos vector [311 [32] and 
the two need never appear in the same basis. 

5.2. Non-Abelian symmetries 

If the symmetry group is large enough that some elements do not commute with each 
other, then it is no longer possible to construct a basis that simultaneously diagonalizes 
all of the generators hence the approach of the previous section needs some modification. 
Instead, we label the basis states by quantum numbers that denote the representation, 
which is no longer simply related to the group elements themselves as the representations 
are in general no longer ordinary numbers, but instead are matrices of dimension > 1, 
and Eq. ( 1241) no longer applies. For SU(2), we choose to label the representations by 
the total spin s, being related to the eigenvalue of of the spin operator, S 2 = s(s + 1). 
Assuming that all of the required operations can be formulated in terms of manipulations 
of these representations, we have a formulation that is manifestly SU(2) invariant; the 
rotational invariance is preserved at all steps and at no time in the calculation is it 
necessary to choose an axis of quantization [36]. This supersedes the earlier approach 
based on the Clebsch-Gordan coefficients [3D]. The non-Abelian formulation is an 
important optimization, because it increases the performance of the algorithm by an 
order of magnitude or more [36] compared with the Abelian case, while enabling more 
accurate and detailed information about the ground-state magnetization. The basic 
ingredient that enables this rotationally invariant construction is the Wigner-Eckart 
theorem [41] , which we can state as: When written in an angular momentum basis, 
each matrix element of an irreducible tensor operator is a product of two factors, a 
purely angular momentum dependent factor (the "Clebsch-Gordan coefficient") and 
a factor that is independent of the projection quantum numbers (the "reduced matrix 
element"). We formulate the algorithm in such a way that we store and manipulate only 
the reduced matrix elements, factorizing out completely the Clebsch-Gordan coefficients. 
The efficiency improvement resulting from the non-Abelian formulation is that the 
matrix dimensions m and M now refer to the number of irreducible representations 
in the basis, which is typically much smaller than the total degree of the representation. 
For a scalar state, this equivalence is precise: a single representation of degree N in 
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the non-Abelian approach results in N degenerate eigenstates when the symmetry is 
not used, with a corresponding improvement in efficiency We do not give here a full 
introduction to the theory of quantum angular momentum, rather we present, in the 
style of a reference, the important formulas required to manipulate MPS wavefunctions 
and operators. For a comprehensive introduction see for example references |l2l H3] . 

Using the normalization convention of Biedenharn [42J, we define the matrix 
elements of a tensor operator transforming as a rank k tensor under SU (2) rotations, 
as 

( j'm' 1 1$ I jm) = ( f || TW I) j ) Cl&C , (27) 

where C.'.'.' is the Clebsch-Gordan (CG) coefficient, j'j label the representation of SU(2), 
and m = —j, + and m' = —j', —j' + 1, . . . , / label the projections of the 

total spin onto the z-axis. Using the orthogonality of the Clebsch-Gordan coefficients, 
this defines the reduced matrix elements, 

( f II TM (I j ) = £ CiXL < j'm' | t£ ] I jm ) , (28) 

mM 

where m! is arbitrary. Note that this normalization is not the same as that used 
by Varshalovich et. al [13], whom instead use an additional factor \/2k + 1 in the 
reduced matrix elements. This is a tradeoff; some formulas simplify slightly with this 
normalization, but the normalization used here has the advantage that the reduced 
matrix elements of scalar operators (with k — 0) coincide with the actual matrix 
elements as all of the relevant Clebsch-Gordan coefficients are equal to unity. Given the 
definition of the reduced matrix elements, we formulate the remaining formulas without 
further reference to the axis of quantization, except as an intermediate step to relate 
the reduced matrix elements prior to factorizing out the Clebsch-Gordan coefficients. 
The coupling of two operators is just as for the coupling of ordinary spins; 

' S [ki] x T M [k] 



(29) 



which denotes the set of operators with components 



c[ki] T [k 2 ]] [k] = v r kik2k ^ ki ^T^ ttrfi 



Applying the Wigner-Eckart theorem gives, after a few lines of algebra, 

(/|| [sNxTN] W (I j) 
= (-1)'+''+* Ej" v / (2/' + l)(2fc + l)| £ k } ] " k | (31) 

x(/||S[ fel l || /'}(/' || 

where {• • •} denotes the Qj coefficient (121 H3j. 

A special case of the coupling law Eq. (13T1) that we will need is when the operators 
act on different spaces, such that they have a tensor product form 

s$p = rw (i)®/(2), 

= /(1)® T W(2). 1 > 
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Here I[i) denotes the identity operator and T^ ki \i) is an irreducible tensor operator 
with respect to the angular momentum J(i) of part i of a two-part physical system 
(z = 1,2). The total angular momentum of the system is J = J(l) + J{2). In this case, 
we write the coupling as [S [kl] xT H ] W = [T [kl] (l)®T [k2] (2)f ] . Repeated application of 
the Wigner-Eckart theorem to these tensor operators gives, after some algebra, 

(f (j'd&ioti II [T^(l)®T^(2)f || j ( 3l32ai a 2 ) > 



3i 32 3 
ki k 2 k 

3i 3 2 3 



[fi K) II T [fcl] (l) || jt (a 1 ))(f 2 K) || rW(2) || j 2 (a 2 ) ) , 



(33) 



where 



3i 32 3 

[(2j[ + l)(2f 2 + l)(2j + l)(2k + 1)]U h k 2 k \, (34) 



3i 32 3 
ji h k 
3i 3 2 3 

and the term in curly brackets is the Wigner 9j coefficient, which can be defined as a 
summation over 6j coefficients [12] 



Ji 32 f 




iyi+j2+j+k 1 +k 2 +k+j' 1 +j! 2 +j'^2^_^2j"^2j" _)_ 1 

j" 

3' h j" \ j f f 2 j[ \ j j" ji f 2 \ 



(35) 



h j k j 1 ji ki f 



32 h j 



We can define an operator norm, corresponding to the usual Frobenius norm, such 



that 



\X [k] || 2 = Tr X [k] ■ X m = Tr X m ■ X [k] 



After some arithmetic, we see that 



l* [fcl HL = £(2j' + i)IO" II t« ||j)| 2 



(36) 
(37) 



j'j 



For the center-matrix formalism, we need the transformation 

A [s] ij ^ E C* A' [s] kj (38) 

k 

where k is a d x m dimensional index that encapsulates both a s' and a j' inde:x[§| : 
k ~ (s',j'). Requiring A'^ s \j to satisfy the right orthogonality constraint, A'A^ = 1, 
this requires 

A' [s] kj = 5 rj 5 s , s [with k ~ (a', /)] (39) 

with 



Cifc 



(40) 



More precisely, k runs over the Clebsch-Gordan expansion of s' <8> j'. 
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In the other direction, we need 

AWy -> J2 A ' [S] ik C V ( 41 ) 

k 

where k ~ (s', i'). Requiring A^ s \k to satisfy the left orthogonality constraint, A'^A = 1, 
this requires 



A' [s \ k = 5 s ,aJ^±± [with k ~ (*', 0] (42) 



and 



The most natural definition for a matrix product operator has two lower indices 
and three upper, 

AfW- (44) 

which transforms as the product of two operators of rank [k], with matrix elements 

( s'q'; j'm' | M] fc l | sg; jm) = ( s'; f || AfW || s; j ) C^, ■ (45) 

Note that the product of an operator and a state requires a contraction of the index s, 
which has the symmetry of over two lower indices, and then shifting the result index s' 
from upper to lower. For 577(2), the required phase factor is (— l) s+fe ~ s ', giving the rule 

BM = (MA) [S ' ] = J2{-l) s+k ~ s 'M [k]S ' s ® A w . (46) 

s 

The action of a matrix-product operator on another matrix product operator is 

X [x] = M [m] N [n] ^ ^ 

which corresponds to the ordinary (contraction) product in the local basis and the 
tensor product in the matrix basis, and therefore results in the product of a Qj and a 
9j coefficient from equations Eq. ( 13T1) and Eq. (1551) respectively. 
For the evaluation of matrix elements, we need the operation 

F /j> = E MS a' S a A 7i B j'j E ij ( 48 ) 

s',s,i,j,a 

On expanding out the reduced matrix elements, we see immediately that the coupling 
coefficient is 

is's .r-.fi* 



*" [a V= E 

a,i,j,k,s,s' 



J S J 

a k a' 



i s' a 



M^A^B^F^j (49) 



Conversely, from the left hand side, 

K a 3 = E E^M'lA^B^ (50) 
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is 



2i' + 1 
2i+T 




E 



a k a' 



(51) 



a' 



s 



I S I' 



■I 



On interchanging A ^ E, B ^ F, this becomes the equation for a direct operator- 
matrix-product multiply. But using the center-matrix formalism, we want instead the 
operation 



j'jk 

where C and C transform as scalars, i.e. the quantum numbers impose i' = i, f = j. 
This is essentially a scalar product E ■ F, and the coupling coefficients drop out. 

6. Conclusions 

In this paper, we have presented an introduction to the MPS formulation of the DMRG 
algorithm for the calculation of ground- and excited states of one- dimensional lattice 
Hamiltonians. The MPS formulation is extremely flexible, allowing the possibility 
for algorithms that act on several distinct wavefunctions at once. The simplest such 
algorithms are for the fidelity and expectation values involving unrelated wavefunctions, 
(■010) and (if)\M\<f>), which are difficult to extract from conventional DMRG. This 
gives access to new tools for the analysis of quantum phase transitions, by measuring 
the scaling function and exponents for the fidelity between ground-states as a function 
of the interaction strength. In addition, the MPS formulation allow optimized versions 
of algorithms for dynamical correlations [2H[32] and time evolution [T7j, which remains 
a fertile for continued algorithmic improvements. 

Finally, we note that in the simulation of finite temperature states via a density 
operator or purification [T2], [H] in the absence of dissipative terms that mix the particle 
numbers between the real and auxiliary systems, the symmetries of the system are 
doubled, such that the symmetries of the Hamiltonian are preserved by the real and 
auxiliary parts independently. For simulations in a canonical ensemble, this leads to a 
significant efficiency improvement that, as far as we know, has not yet been taken into 
consideration. 
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